library(sasld)

# -------------------------------------------------------------------
# ATTENZIONE: la funzione sample e' stata modificata 
# a partire da R 3.6.0 per cui i "vecchi risultati"
# non posso essere riprodotti di default
# http://hostingwin.unitn.it/micciolo/PeM/01-rng.pdf
# -------------------------------------------------------------------
# se si vogliono riprodurre i risultati dei Laboratori
# con le versioni pił recenti di R va eseguita la seguente istruzione
# -------------------------------------------------------------------
RNGkind(sample.kind="Rounding")
# -------------------------------------------------------------------


# -------------------------------------------------------------------
# Sezione 15.1 - Il test di Wilcoxon per due campioni indipendenti
# -------------------------------------------------------------------
job <- function(x) setdiff(c(1:5),x)
X <- t(combn(5,3))
Y <- t(apply(X,1,job))

X[3,]
Y[3,]

sA <- apply(X,1,sum); sB <- apply(Y,1,sum) 
mA <- apply(X,1,mean); mB <- apply(Y,1,mean) 
delta <- mA - mB
tmp <- cbind(sA,sB,mA,mB,delta)
colnames(tmp) <- c("Somma A","Somma B","Media A","Media B","Diff.")
tmp


plot(prop.table(table(round(delta,2))))


choose(64,32)


#
lozione <- c(2,4,5)
uv <- c(1,3)
wilcox.test(lozione,uv,exact=TRUE)
wilcox.test(uv,lozione,exact=TRUE)


# Il test di Wilcoxon "approssimato"
lozione <- c(2,4,5)
uv <- c(1,3)
wilcox.test(lozione,uv,exact=FALSE,correct=FALSE)

wilcox.test(lozione,uv,exact=FALSE,correct=TRUE)


#
data(rt)
rnk <- rank(rt$rt)
aa <- tapply(rnk,rt$gruppo,mean)
bb <- tapply(rnk,rt$gruppo,sum)
cc <- cbind(aa,bb)
rownames(cc) <- c("Cellulare","Controllo")
colnames(cc) <- c("Media dei ranghi","Somma dei ranghi")
cc

cellulare <- rt$rt[rt$gruppo=="cell_phone"]
controllo <- rt$rt[rt$gruppo=="control"]
wilcox.test(cellulare,controllo,exact=FALSE,correct=FALSE)


# Il test di Wilcoxon "quasi esatto"
cellulare <- rt$rt[rt$gruppo=="cell_phone"]
controllo <- rt$rt[rt$gruppo=="control"]
y <- c(cellulare,controllo)
x <- c(rep(1,length(cellulare)),rep(2,length(controllo)))
rnk <- rank(y)
tapply(rnk,x,sum)

set.seed(123456)
yy <- sample(rnk,64,replace=FALSE)
tapply(yy,x,sum)

#
nrep <- 100000
out <- numeric(nrep)
set.seed(123456)
for (i in 1:nrep) {
    yy <- sample(rnk,64,replace=FALSE)
    out[i] <- tapply(yy,x,sum)[1]
}

ok <- which(out >=1215); length(ok)
length(ok)/nrep*2

wilcox.test(cellulare,controllo,exact=FALSE,correct=FALSE)


# Un test non parametrico di permutazione 
# e un intervallo di confidenza bootstrap
y <- c(cellulare,controllo)
x <- c(rep(1,length(cellulare)),rep(2,length(controllo)))
ymn <- tapply(y,x,mean)
c(ymn,ymn[1]-ymn[2])

#
set.seed(123456)
yy <- sample(y,64,replace=FALSE)
ymn <- tapply(yy,x,mean)
c(ymn,ymn[1]-ymn[2])

#
nrep <- 100000
out <- numeric(nrep)
set.seed(123456)
for (i in 1:nrep)
{
    yy <- sample(y,64,replace=FALSE)
    ymn <- tapply(yy,x,mean)
    out[i] <- ymn[1] - ymn[2]
}
mean(out)


#
ymn <- tapply(y,x,mean)
(dlt <- ymn[1]-ymn[2])

ok <- which(abs(out) >= dlt)
length(ok)/nrep


# bootstrap
nrep <- 100000
set.seed(123456)
ycel <- sample(cellulare,nrep*length(cellulare),replace=TRUE)
ycon <- sample(controllo,nrep*length(controllo),replace=TRUE)
ycel <- matrix(ycel,nrow=nrep)
ycon <- matrix(ycon,nrow=nrep)
ysta <- apply(ycel,1,mean) - apply(ycon,1,mean)

quantile(ysta,prob=c(0.025,0.975),type=1)


# Differenza fra due mediane
cellulare <- rt$rt[rt$gruppo=="cell_phone"]
controllo <- rt$rt[rt$gruppo=="control"]
a <- data.frame(n=c(length(cellulare),length(controllo)),mediana=c(median(cellulare),median(controllo)))
rownames(a) <- c("Cellulare","Controllo")
a

wilcox.test(cellulare,controllo,conf.int=TRUE,exact=FALSE,correct=FALSE)


#
D <- outer(cellulare,controllo,"-")
D <- as.vector(D)
median(D)




# -------------------------------------------------------------------
# Sezione 15.2 - Il test di Kruskal-Wallis
# -------------------------------------------------------------------
rare <- c(1.75, 3.15, 3.50, 3.68)
occasional <- c(2.00, 3.20, 3.44, 3.50, 3.60, 3.71, 3.80)
regular <- c(2.40, 2.95, 3.40, 3.67, 3.70, 4.00)
y <- c(rare,occasional,regular)
x <- c(rep("rare",4,),rep("occasional",7),rep("regular",6))
rnk <- rank(y)
a <- data.frame(n = c(length(occasional), length(rare), length(regular)), mediana =
                    c(median(occasional), median(rare), median(regular)), rango_medio = tapply(rnk,x,mean))
rownames(a) <- c("occasionalmente","raramente","regolarmente")

kruskal.test(y ~ as.factor(x))




# -------------------------------------------------------------------
# Sezione 15.3 - Il test dei segni
# -------------------------------------------------------------------
data(georgia)

internet <- georgia[,11]
tv <- georgia[,12]
delta <- internet - tv
table(sign(delta))


prop.test(19,(19+35),p=0.5,correct=FALSE) 

binom.test(19,(19+35),p=0.5) 

prop.test(19,(19+35),p=0.5,correct=TRUE) 




# -------------------------------------------------------------------
# Sezione 15.4 - Il test di Wilcoxon per due campioni dipendenti
# -------------------------------------------------------------------
prima <- c(2.5,4.0,1.5,4.5,2.0,4,3,3.5,3.0,2.0,1.0,4)
dopo  <- c(3.0,3.5,3.0,5.0,4.5,5,5,2.5,3.5,4.5,2.5,5)
delta <- dopo - prima
table(sign(delta))

binom.test(10,12)

rnkd <- rank(abs(delta))
rbind(delta,rnkd)

table(delta)

set.seed(123456)
(sgn <- sample(c(-1,1),12,replace=TRUE))
(y <- sgn * rnkd)

ok <- which(y > 0); sum(y[ok])


#
nrep <- 100000
out <- numeric(nrep)
set.seed(123456)
for (i in 1:nrep)
{
    sgn <- sample(c(-1,1),12,replace=TRUE)
    y <- sgn * rnkd
    ok <- which(y > 0)
    out[i] <- sum(y[ok])
}


y <- rnkd*sign(delta)
ok <- which(y > 0); sum(y[ok])

ok <- which(out >=69.5); length(ok)

length(ok)/nrep*2

wilcox.test(dopo,prima,paired=TRUE,exact=FALSE,correct=FALSE)




# -------------------------------------------------------------------
# Sezione 15.5 - Un confronto fra la potenza del test t e quella del test di Wilcoxon
# -------------------------------------------------------------------
power.t.test(delta = 2, sd = 2.1, sig.level = 0.05,
             power = 0.8, type = "two.sample",
             alternative = "two.sided", strict = TRUE)

power.t.test(n = 19, delta = 2, sd = 2.1, sig.level = 0.05,
             type = "two.sample",
             alternative = "two.sided", strict = TRUE)


#
set.seed(123456)
a <- rnorm(19,mean=80,sd=2.1)
b <- rnorm(19,mean=82,sd=2.1)
t.test(a,b,var.equal=TRUE)

wilcox.test(a,b)


#
nrep <- 100000
n <- 19
ds <- 2.1
set.seed(123456)
aa <- rnorm(19*nrep,mean=80,sd=2.1)
bb <- rnorm(19*nrep,mean=82,sd=2.1)
A <- matrix(aa,ncol=n)
B <- matrix(bb,ncol=n)

ris <- matrix(nrow=nrep,ncol=5)
for (i in 1:nrep)
{
    a <- A[i,]
    b <- B[i,]
    vara <- var(a)
    varb <- var(b)
    varp <- (vara + varb)/2
    tt <- t.test(a,b,var.equal=TRUE)$p.value
    wt <- wilcox.test(a,b)$p.value
    ris[i,] <- c(mean(a),mean(b),varp,tt,wt)
}

colMeans(ris[,1:3])

table(ris[,4] < 0.05)
table(ris[,5] < 0.05)
table(ris[,4] < 0.05,ris[,5] < 0.05)